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Abstract. In late sixties, Mihail Budyko and William Sellers, a Russian and an American 
climate scientists, independently introduced the concept of Energy Balance Model with ice albedo 
feedback. Since then many have followed in their footsteps to establish various versions of this 
model. In this paper, a novel equation is introduced to account for the dynamics of the ice line, and 
is coupled to Budyko's model. We found that the coupled temperature profile-ice line system has 
a one dimensional center stable manifold. Furthermore, this conceptual model, that accounts only 
for feedback from the ice albedo and that should not be used for predictive purposes, suggests that 
under some climate forcings, the ice free earth is unstable. 
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1. Introduction. Energy Balance Models (EBMs) rest on the concept that in 
the earth's equilibrium state the energy received from the sun's radiation is balanced 
by the energy re-emittcd back to space at the earth's temperature. Budyko's EBM 
aims to model the latitudinal distribution of the surface annual mean temperature, 
by taking into account the feedback from ice albedo (icecap's reflectivity factor). As 
a dynamical system then, Budyko's model takes place in an infinite dimensional state 
space. While the discussion of climate feedbacks easily enters realm of great complex- 
ity, the concept of ice albedo feedback and its effect on the planet, is not difficult to 
explain. Ice albedo is the degree of which ice reflects light or energy. Since ice cover 
in general has whiter color than the blue ocean, it reflects more light or energy back 
into space. This creates a positive feedback, properly called the ice albedo feedback, 
because more ice means less energy absorbed, which induces cooler climate favoring 
condition for more ice formation, and so on. On the other hand, shrinkage of ice 
mass induces more ocean surface, which absorbs more energy, creates warmer climate 
and in turn, promotes more melting. But despite this easy to explain concept, the 
conclusions on the stability of icecaps vary greatly, even among the simplest mod- 
els. While previous results were useful in understanding steady states of the energy 
balance models, they did not include mechanisms for ice-water boundary or iceline 
movement, making it difficult to obtain rigorous results for icecap stability. In this 
paper, to go along with Budyko's model, we introduce an equation that describes the 
movement of the ice line. We call this coupled system. Dynamic Iceline Budyko's 
Model (DIBM), and it consists of two equations: the first is a version of Budyko's 
model, similar to that written by Ka Kit Tung [8] describing the evolution of the tem- 
perature profile, and the second is the novel equation that models the iceline dynamics. 

Our planet currently has a small ice cover surrounding both poles, and in light 
of the ice albedo feedback concept, many has pondered on the stability of these po- 
lar ice caps. The idea of small ice cap instability, or sometimes abbreviated SICI, 
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perhaps was started in 1924 by CEP Brooks [14], [13], when he hypothesized that 
in this planet, because of a mechanism later coined as ice albedo feedback, only two 
stable climates are possible: ice free and large ice capped earth, with the cryospheres 
extending from the poles pass the 78° latitude. In late sixties, a Russian climatologist 
Mihail Budyko proposed an ice albedo feedback model based on the energy balance 
principle, where transport is represented by a simple relaxation process in which the 
temperature of each latitude dissipates to the global average temperature. He con- 
cluded that both ice free and polar or small ice cap climates are unstable (see p. 618 
[7]). Around the same time, William Sellers from University of Arizona proposed a 
minimal complexity climate model which includes planetary albedo, but uses separate 
atmospheric and oceanic transport processes [12]. In the next decade, Gerald North 
explored in depth the diffusive transport version of the energy balance model with 
the ice albedo feedback, and in his works, upon varying some of the parameters, the 
small ice cap instability disappears [2], [3]. 



While the small ice cap instability discussions in the mid twentieth century were 
fuelled by the possibility that the planet is going into its glaciation period, the more 
recent discussion on the state of the cryosphere is motivated by the rapid decrease of 
the Arctic sea ice extent, induced by global warming and placed more emphasis on 
the reversibility of an ice free state in a warmer climate. Recent results that attempt 
to answer this question have also been characterized by the use of computer simula- 
tions such as the one developed by in 2008 by Merryfield, Holland and Monahan [TT] . 
a simple climate model based on the Community Climate System Model version 3 
(CCSM 3)-a global climate model developed by NCAR. The simple model is nonlin- 
ear and admits abrupt sea ice transitions resembling those in the CCSM 3 simulations. 
In early 2009, Eisenman and Wettlaufer 1 examined an energy balance model with 
seasonal features and a nonlinear forcing given by the sea ice thermodynamics. Their 
results suggest that there exists a critical threshold in the climate warming, beyond 
which a sudden loss of the remaining seasonal sea ice is possible. Another recent 
paper published in 2009 by Notz [13| summarized the discussions on the state of the 
" existence of the cryospheric tipping points" . 



We will show that Dynamic Iceline Budyko's Model admits an unstable large ice 
cap, a stable ice covered planet, a stable small ice cap, and an unstable ice free planet 
indicating polar ice cap loss reversibility. Furthermore, we show that the infinite di- 
mensional system has a 1-D center stable manifold, hence, the dynamic is essentially 
one dimensional. The existence of this one dimensional attractor will also be illus- 
trated in the computer simulations, but in the end, it is the mathematical analysis of 
the dynamics that gives us the confidence in the numerical executions. 



Unlike the more complex computer simulation models, our result should not be 
seen as predictive, but instead, as a tool to understand qualitatively the mechanism 
involved in climate processes. The advantage of our formulation is that its tractability 
allows us to isolate the effects of a single climate feedback, ie, the ice albedo feedback. 
For example, we see from this result that despite being a positive feedback, the ice 
albedo feedback is not enough to cause a tipping point phenomena as seen in recent 
studies, see [T] and reference therein. The comparison of these two results indicates 
the fruitfulness of exchanges between computer simulation models and theoretical 
analysis, hence, between climate scientists and mathematicians. 
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We will organize this paper as follows: in section 2, we introduce Budyko's energy 
balance model along with a new feature to account for the iceline dynamics, and we 
discuss the equilibria of the systems. Then we show some animations that results 
from the numerics to illustrate the dynamics. Section 3 discusses the analysis of 
the Dynamics Iceline Budyko's model and its map, with detailed proofs using graph 
transform method in Section 4. Section 5 provides concluding remarks and a sketch 
of some future directions. 

2. Dynamic Iceline Budyko's Model. The equation governing the temper- 
ature profile evolutions by itself does not induce the movement of the iceline, as 
illustrated by computer simulations in the following section. To capture the feed- 
back from the ice albedo, we propose a system of two time scales, integro-difference 
equations governing the temperature distribution and the iceline dynamics. With this 
proposed system, we asked the following questions: 

• What is the appropriate function space for the temperature profiles? 

• What are the dynamics of the proposed system? 

• What are the parametric conditions to have an equilibrium state at the cur- 
rent ice line? 

• What does the model suggest if the current condition is perturbed so that 
the planet is ice free? 

The main results of this study are: 

• The identification of the function space in which the proposed system is well 
defined. 

• The existence of an invariant manifold for the Euler approximation of the 
proposed system. 

• A parametric condition to have an equilibrium at the present climate. 

• The instability of the ice free earth. 

We first introduce the details and some background of the model, then we explore 
the equilibria, the invariant manifold and, finally, the instability of the ice free planet. 

2.1. Details of Budyko's Model. While Budyko introduced his original model 
as a concept, following his steps many have formulated and parametrized the energy 
balance and ice albedo feedback concept. KK Tung has summarized the formulation 
as a differential equation with a nice exposition in his book Tung (2007): 

R^T{y) = Q ■ s{y) ■ [1 ~ a{y)] ^[A + BT{y)] - C ■ [T{y) - J^' T{y)dy] (2.1) 

The function T[y) is the annual average surface temperature at the zone y, and 
its graph over the domain of y is called the temperature profile. Other previous 
authors, eg. [2], |10j . also treated the model as a differential equation, hence, presup- 
posed the continuous dependence of the annual average temperature profile T{y){t) on 
the time variable t. We argue that the evolution of the yearly averaged temperature 
profile T should be governed by a difference equation. 
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Let At[/](fc) denotes the difference equation of / at time node k, tliat is 



We consider tlie following Budyko's difference equation at each time node k: 
RAt[Tiy)]{k) = Qs{y){l - a{Tj){y)) - {B + C)T{y){k) - A + CT{k) (2.2) 



Here, the constant R measures the heat capacity of the surface, and we assume 
that i? = 1. Such assumption does not change the qualitative behavior of the system. 
The independent variable y is sin(0) , where 6 represents some latitude in the northern 
hemisphere, therefore y lies on the unit interval. This model assumes that the annual 
distribution of the surface temperature is symmetric about the equator. 

We will now examine in details the terms on the right hand side of the EBM 
equation above. The constant Q is the solar constant, representing the amount of 
energy received from the sun at the top of the atmosphere. The function s{y) is the 
latitudinal distribution of that energy which could be computed from earths orbital 
elements. Many authors, such as KK Tung [S], and North [2] used the Legendre poly- 
nomial function approximation s{y) = 1 — 0.482 '^^^""'^ for this distribution. We will 
use the same approximation as well. 

The term 1 — Oi{y) represents the fraction of the radiative energy absorbed by the 
earth at location y. To establish the ice albedo feedback effect, it is assumed that ice is 
formed when the temperature at a certain location stays below a critical temperature 
Tc, taken to be — 10°C. Also, it is assumed that the surface is either water (ice free) or 
ice covered and that there is only one ice line, 7]. Since ice reflects more sunlight than 
water does, area covered with ice has higher albedo than that covered with water. In 
this paper we use the approach that the albedo function a{y) at location y depends 
on y relative to the location of the iceline, and not on the temperature T{y). Let 
77 G [0, 1] denote the location of the ice line. The albedo function we use in this paper 
is smooth and is iceline dependent: 

Definition 2.1. Iceline- Dependent Smooth Albedo 

Given that the ice line is at rj, the albedo at y is 

a{v)iy) = 0.47 + 0.15 • tanh[M ■ {y - 77)] 

Here M is the parameter representing the steepness of the albedo near the icelin- 
eand is a fixed quantity, presumably dictated by the planet. 

As in Tung 8J , balancing the absorbed radiative energy contained in the first term 
are the re-emission term, A + BT{y) and the transport term, C ■ {T{y) — T), where T 
is the global average J,„ T{y)dy. The constants A = 202 watts per squared meters 
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and B — 1.9 watts per squared meters per °C are derived from fitting a linear func- 
tion through satehte data of Outgoing Long-wave Radiation (OLR) at the top of the 
atmosphere [5^. The constant C in the transport term is taken to be 1.6B = 3.04 and 
is chosen so that one equihbrium fits the current chmate with ice hue near the pole. 



2.2. Dynamic Iceline. So far we have an equation that describes the evolution 
of the temperature profile over the northern hemisphere, and that takes into account 
the ice albedo feedback. We will add to this an equation that describes the evolution 
of the ice line 77. Previous literatures such as Budyko [7], Tung [Rj, North [2], etc. have 
used the idea that ice is formed when the temperature is below a critical threshold, Tc- 
We adapt this idea to the ice line evolution by prescribing a poleward movement of 
the ice line when the temperature at the ice line T{ri) is above a critical temperature 
Tc, and equatorward movement otherwise. Also, we assume that the movement of the 
ice line happens at a much slower rate compared to the evolution of the temperature 
profile. Therefore, for e << 1 and with Tc = —10 as used by previous authors, the 
equation for the ice line evolution at time node k can be written as: 



Advm = "^^ + ^1 ^^^^ = . • [Tivm - TJ (2.3) 

Some recent computation has suggested that to be physically relevant, the value 
of e is very small, possibly in the order of 10^^ [22]. We will see in the simulation 
that using a much larger e = 0.025 the attracting 1-D manifold still appears. 

The object of our analysis, Dynamic Iceline Budyko's Model, is therefore the fol- 
lowing infinite dimensional, two time scale, integro-difference equations: 



|A4T(y)](fc) = T([T(y),r7])(fc) 

with F and G as the following: 

F{[T{y),r,]) = Qs{y){l ~ a{r,){y)) - {B + C)T{y) -A + C T T{y)dy 

Jo 

Gmv)M)^<nri)-Tc) 

2.3. Equilibria and Animations. When the ice line 77 is fixed, the tempera- 
ture profile equilibrium, with ice line at r] is: 

w , g ■ sjy) ■ (1 - a{rj){y) - A + C T* {rj){y)dy 

T my) = 

We call the set of the Lipschitz continuous functions T* := {T*{rj) : 77 G [0,1]} the 



local equilibrium set. This set is the steady state of the first equations in (2.4). 
Notice that, without a mechanism specifying the movement of the iceline, the tem- 
perature profile will dissipate to the local temperature profile T*{ri){y) at the initial 
ice line rj. The following figures simulate the evoluti ons of the temperature profile 
To(j/) 14 — 54y^ following only the first equation of (2.4 1, with the ice lines fixed at 
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% = 0.1,0.3,1. 

Note: on the electronic copy, the following figures are the initial temperature 
profile and the initial ice line, and the local equilibria at the respective ice lines. Click 
anywhere on the left figures to start the animations following Budyko's equations. 
Similar animations can also be found on the web through: |http : //math . arizona . 
,edu/~ewidiasih/ index . html/ewidiasih/Research . htmlj 




Fig. 2.1: The horizontal axes are the domain of the temperature profile T. The solid 
blue curve is the temperature profile, and the big red 'X' represents the location of the 
ice line. The figures on the left are the initial temperature profile, T{y) = 14 — 54?/^ 
and ice lines at ?7 = 0.1,0.5 and 1, and the figure on the right is the steady state 
temperature profile T*(r/)(y) at the respective ice lines. Clicking on the top figures 
will start the animation. 
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The next animations track the dynamics using the proposed system with param- 
eter e = 0.025 and with the same starting temperature profile and ice hues as the 
previous series of figures. It should be noted here that we choose this value of e, not 
based on the physical considerations, but rather to create a reasonable animations run. 

Notice that the initial temperature profile first evolves to a temperature profile 
similar to the local equilibrium temperature profile, one with a large slope at some 
ice line 77 close to the initial ice line, and eventually they together move toward an 
equilibrium. This suggests the existence of an invariant manifold in the phase space 
(T, 77) . We can directly compute the global equilibria for the proposed system by 
setting the ice line temperature of the local equilibria to the critical temperature and 
they are: 

(r*(??i)(y),??i) where 771 ^ 0.225, and (T*{r]){y),r]2) where 772 = 0.962 




Fig. 2.2: The figures on the left are the initial temperature profile, T(y) = 14 — 54?/^ 
with ice lines at 770 = 0.1, 0.5 and 1, and the figure on the right is the steady state 
temperature profile T*(?7)(j/), with the ice lines at the equator or ij ~ 0, rj ^ 0.962. 
As before, clicking on the top figures will start the animations. 
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3. Analysis of Dynamic Iceline Budyko's Model . Our goal is to explain 
the dynamics of the temperature profiles T{y), coupled with the ice line rj in the 
simulation above. The main result is the existence of a one dimensional center stable 
manifold. The idea of an attracting invariant manifold in a fast slow system has been 
explored and developed by many including J. Carr 18J, S. Wiggins [17], Fenichel [20] . 
and Vanderbauwhede |19j to name a few. 

First, we argue that it is reasonable to assume that the albedo functions has 
bounded local variation, that is, for each fixed ice line r], the albedo function a{ri) has 
a bounded Lipschitz contant. Let M be this bound, ie. for any real values x and z, 

|a(77)(a;) - a{'q){z)\ < M\x - z\. 

When the ice line is at the equator, we have an ice covered planet, otherwise, 
when it is at the pole, our planet is ice free. The forcings in Budyko's model as 
explained in section 2 is fixed to the current climate which allow for a small ice cap 
as an equilibrium. As mentioned in the introduction, an important question in light 
of the current global climate issue is: 

"What does the model suggest if this equilibrium is perturbed so that the planet 
suddenly becomes ice free?" 

To be able to describe the dynamics of the ice line ry at the end points, we will 
embed the ice line interval and the domain of the temperature profile in the real line. 
We will also embed the vector field in such a way that preserves the dynamics in the 
unit interval including its end points. This embedding allows the dynamical analysis 
of the polar and the equatorial ice line and is a mathematical convenience for showing 
the existence of an inertial manifold. 



We thus consider the following system similar to (2.4 1 



MT{y)] 



G([r(y),77]) 



(3.1) 



where we define: 
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BT{y)] + C[J^ T{y)dy - T(y)], when y < 
BTiy)] + C[J^ T{y)dy - T(y)], when < y < 1 
BT{y)\ + C{\1 T{y)dy - T(y)], when 1 < y 



and 



G([T(y),7,]):=e[r(,7)-T,] 



Notice here that for a small e, the vector field H separates the dynamics into 2 
time scales: the fast dynamics for the temperature profiles, and the slow time scale 
for the ice line evolution. 



Also, observe that in the extended version, any zone y outside of the unit interval 
has the same solar forcing as the closest endpoint. The equilibria of the extended 
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version are similar to the original version. The local equilibria for the temperature 
profile with ice line at 77 is the following Lipschitz continuous function: 



Q.s(O).(l-a(.0(O»-A+Cj„^T-(,)fa) ^^^^ ^ ^ ^ 

T*iv)iy) = { Q-siyyi^-<nKy))-A+cf„^T^n(y)) ^^^^ 0<y<i 
Q.siiyii-ai,mj-A+cf„^T'inKy) ^^^^^^ ^ ^ ^ 

3.1. Results. The interval of the ice line is M, and we take the function space 
for the temperature profiles T{y) to be the Banach space 

i3 = {T : M M : r is bounded and continous, with bounded Lipschitz constant} 

with the norm ||r||e = ||T|U 

3.1.1. Inertial Manifold. We now consider the following shift operator to Dy- 



namic Iceline Budyko's Model (3.1) in B 
mt([r,77](fc)) := 



T{k + t)^Tik)+t-F{[T,T^]{k)) 
rjik)^rj{k)+t-Gi[T,rj]{k)) 



We obtained the following results: 

Theorem 3.1. For an e small, there exists an attracting invariant manifold for 
the shift operator of Budyko 's equation that is, 



1. There exists a Lipschitz continuous map 



2. There exists a closed set V d B, such that for any {Tq, 770) G 2? x M, the distance 
dist[mt{[TQ,r]Q]){k),{^*{ri),ri)] decreases exponentially as k increases. 



Let $* denote this invariant manifold. As a corrolary to the theorem (3.11, for 
r] £ [0,1] we can compute the distance of the invariant manifold $*(?7) to the local 
equilibria manifold T* and more importantly, we can describe the asymptotic dynam- 
ics of the ice line. The following graph is the graph of ice line temperature of the local 
equilibrium in the extended version, T*{vi)[rf) over the interval [0, 1]. 

Many authors uses the critical temperature T^ = —10. The dynamics of the ice 
line is determined by the sign of the difference between the temperature at the ice 
line and the critical temperature T{r]) — Tc, which is a one dimensional function of the 
ice line rj. To determine the dynamics of the ice line at the steady state temperature 
profile, we graph T*{ri){ri) - T^ over the interval [0, 1]. The function T*{r]){r]) - 
has two roots in the unit interval, as in the previous section, let the root closest to 
the equator, ?; = be 771 and the other root be 772. 
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We conclude from Theo rem |3.1| that the dynamics reduces to 1 dhnension. This 
fact is ihustrated in Figure (3. lb I. Here, we observe that when the temperature pro- 



file reaches a steady state, any ice line rj £ (?yi,?72) moves to the right because the 
ice line temperature T*{ri){ri) exceeds the critical temperature. On the other hand, 
when < ?7i, and t] > r]2, the ice line temperature T*{ri)[ri) falls below the critical 
temperature, and therefore the ice line rj moves to the left. In particular, we note 
that when the ice line is at the pole, that is when 77 = 1, then this ice line advances 
toward the equator. 



Corollary 3.2. In the unit interval, the invariant manifold $* is within 0(e) 
of the local equilibrium set T* . Therefore, for a small e, the ice free planet is unstable. 

As an example, for a small e, the equilibrium iceline temperature ^*{r]){ri) is 
within 1°C of T*{rj){ri), the ice free earth is unstable. 




n(t) 



Eq. 



BIG ICE CAP 
UNSTABLE 



SMALL ICE CAP 
STABLE 



Fig. 3.1: Figure on the left illustrates <&*(?/)(??) within \°C of T* (77) (77). Figure on the 
right is the reduced 1-D dynamics 



In light of the current global climate discussion, this result suggests an optimistic 
view. Suppose that the planet's temperature somehow rises quickly that the ice line 
moves to the pole and the small ice cap disappear. If the climate parameters, ie. 
the parameters A, B, C are kept the same, then the temperature profile will dissipate 
toward that on the invariant manifold, so that the ice line temperature falls below the 
critical temperature, and the ice line will start moving toward the stable equilibrium, 
the small ice cap. 

We obtain the following diagram, which is a bifurcation diagram of the solar ra- 
diation parameter Q against the ice line 77 in the equilibrium state. The dash blue 
line denotes unstable regime, and the solid line is stable, while the vertical green line 
is the current solar radiation used, Q = 343. We see that the ice free state ie. ice line 
at the pole ov rj = 1, belongs to the unstable regime when the solar radiation Q is 
near the current value, 343 watts per meter squared. 
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Fig. 3.2: Bifurcation diagram for Q-thc solar input 



4. Technical Details. This section present the proof of the existence of the 1- 
D center stable manifold. Readers with little interest in the mathematical treatment 
may wish to skip this section and proceed to the discussions on some future directions 
and conclusion in the next section. 



We consider the extended version of the Dynamic Iceline Budyko's EBM (3.1 1 



4.1. A Function Space for Temperature Profiles. We take the space for 
the the temperature profiles to be the Banach space: 



:= {T : M — > M : r is bounded, continuous with bounded Lipschits constant} 

with the norm sup norm ||T||oo 

4.2. Equilibria. Lemma 4.1. Two Equilibria 
The system of the differential equations (3.1) has two equilibria. 

Proof. First, we fix 77 and set F{T^ rj) to zero and solve for T, in which T = T{ri){y), 
a function depending on 77 and y. Let T*{r/){y) denote this solution. 



Q..(0).(l-a(,)(OWJ„^T-(,)fa) ^^^^ y ^ Q 

T^mv) - { Q-^(y)<^-^in)iy)-^+cf,]T'My)) ^^^^ o<y<i 

Q..(l).(l^a(,)(l)-A+Cj„^T»(,)fa) ^^^^ ^ ^ ^ 

Let 5(77) := Q ■ s(y)(l - a{y,rj))dy, then T*{rj){y)dy = and substi- 

tuting this to T*(r/) we get an expression for T*(rj) which only depends on rj and y. 
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T*{v){y) = 



Q.s(0).(l-a(,))(0))-^+g(g(r))-A) 

B+C 

Q.s{yy{l-aiv){v))~A+§{g(n)-A) 
B+C 

Q.s{iy{l-a{ri){l)-A+^{gir,)-A) 



, when y < 
, when < y < 1 
when 1 < y 



B+C 



Then we use the second equation to find the ice edge(s) equiUbria, that is, we set 

T*{ri){rj) — Tc to zero, and solve for rj. Let h(ri) denote the function T*{rj){ri) — Tc- We 
can eliminate the case that 77 < or 77 > 1, since the temperature profiles for these 
ice lines do not intersect the critical temperature = —10. 

As mentioned in the introduction, we also assume that the parameter of the 
albedo function's maximum slope, ie. the parameter M is large, eg. M > 10. And 
we need the following computational lemma: 

Lemma 4.2. For M > 10, T*{rf){ri) — has two roots on the unit interval. 

Let rji and 772 be these roots, and suppose that 771 < ri2- The equilibria for the 
above differential equations are therefore {T*,r]\) and (T2,772), with 



4.3. Inertial Manifold. By using Hadamard's graph transform method we will 
show that in an appropriate function space, for e small. Dynamic Iccline Budyko's 
Model has an inertial manifold. The main theorem that we will prove in this section 
is the following: 

Theorem 4.3. Existence of Local Inertial Manifold 
For e small, there exists a locally attracting invariant manifold for the shift operator 
of Budyko 's Energy Balance Model. 

4.4. 2. The space of graphs.. For each ice line 77 G R, we consider (77, 
where we take the space of the graphs $ as the space G of the bounded continuous 
functions, BC^{R,B) with norm \\^\\g = supr,eR\\^{'n)B\\- We will need the following 
lemma: 

Lemma 4.4. 

The set Q with the norm \ \- \ \g is a Banach space. 

Lemma 4.5. 
For each fixed L, the set 



n{y) := T*{m){y) 



n{v) := T*{ni){y) 



Q ■ (1 - a{m){y)) ■ s{y) + 1.6g(77i) - 2.6A 
B + C 

Q ■ (1 - a{y, m){y)) ■ s{y) + l-6g(??2) - 2.6^ 

B + C 



□ 



gL = {^&g: VC,r? e K, mo - $(77)1 Is < L\Q - r?]} 



is a closed set in the Banach space Q. 
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4.5. The albedo function 0(77) (y) as a graph. Given the location of the ice 

line 77, the albedo function of the planet is the fimction a{r])(y) — A7 + .15-tanh{M ■ {y — 
r])). The parameter M is the steepest slope of the function, which occur on the ice line 
T]. Furthermore, the albedo function a as a function of ry is , with the upper bound 
for the first derivative being .15M, therefore, supr,xeK.S''J'Pye«- < .15M. 

Also, for any rj and y, 0.32 < a(r/)(y) < 0.62. Therefore, if L > max{.62, .15M}, 
then the albedo function a is a member oiGL- We will use this bound as the Lipschitz 
bound of the space of graphs. To be precise, we define: 

Definition 4.6. Lipschitz bound, L 
The Lipschitz Bound for the space of graphs is a number L such that 

r > « Hv)-<0\ 
L > swp„,ceR — 

In particular, L > max{.62, .15m} . 

Definition 4.7. Temperature profile bound, r 
The bound for the space of temperatures is the number r such that 

r > supy^io,i]\Q ■ s{y)\ 



We now consider the the action of Dynamics Iceline Budyko's Model over the set 
Gl n B{0, r) in the space of graphs, where B(0, r) = {$ e : 1 1$| | < r}. 

4.6. The Shift Operator of Budyko's Model as a Graph Transform. We 

will define m a transformation that extends the vector field {F, G) so that its action 
is defined for all ice boundaries in the real line and for all temperature profiles in B. 
But, first, we need a lemma: 

Lemma 4.8. If e < j^, then for any fixed t < 1, and for each e M and each 
^ &Ql such that W'^Wg < r, there exists a ^ e M such that 

V = ^ + e-t-mm-Tc) 

Proof. We will show that given $ G ^Jl, there exists a unique k = kis> G BC^{M.) 
such that for any given r], r] = {ri + k{ri)) + e($(r? + k{ri)){r] + k{r})) — Tc). 

We define a transformation T on SC"^(M) by 

{Tk){ri) = et[Tc - $(?? + k{'n)){7j + ^(77))] 

Clearly, | \Tk\ \ < 00. Wo will now show that T is a contraction mapping on the Banach 
space BC''\M.). Indeed, for any two bounded continuous functions ki and k2 with sup 
norm less than r, we have that: 
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\\{Tk,){rj)-iTk,)m 

= et\<i>{i] + ki{r])){T] + ki{ij)) - $(77 + k2{v)){v + k2{v))\ (4-1) 

< ei[|$(ry + fci(?7))(r7 + hiv)) - + k2{r]))\{v + ki{r,)) (4.2) 
+ Mv + k2{v)){v + hiv)) - ^{V + k2{v))iv + k2{v))\] (4.3) 

< etlimfi + kM) - + k2{ri))\\B + m^l + fc2(r/))|b • 1^1(77) - k2{ri)W (4.4) 

< et{L + r)\k^{'n) - fc2(?7)| < p\ki{'n) - k2{ri)\. (4.5) 
with the number p = et{L + r) < 1. 

Therefore, the transformation T is a contraction on a Banach space, and by the 
Banach Fixed Point Theorem, there exists a unique fixed point. Let k be the fixed 
point, that is, 

fc(r/) = k^{Tj) = e(Te - $(7/ + k{7^))). 

Therefore, 

77 + k{r]) =7] + €t{Tc - + k{ri)){ri + k{ri))). 
Letting = rj + k{rj) finishes the proof. □ 

Definition 4.9. Graph Transform using Budyko's EBM 
Given a graph $ G we define Budyko Graph Transform m as: 



i($) =m($)(?7) = $(0+i-i^(«',C) 



with ^ as in Lemma { 4 



To prove the existence of the inertial manifold for Dynamic Icehne Budyko's 
Model, we first show that m is a contraction on Qj^ n 5(0, r). 

Lemma 4.10. As denoted above, let L be the Lipschitz bounds, r be the tempera- 
ture profiles bounds, and B be the the constant from the re- emission term of DIBM. 
There exists an e — e{L, r, B) > small, such that for any fixed time step t, < t < 
, the map m is a contraction on Ql<^ B{0, r). 

Proof. 

We will show that for any fixed time step t < -g^, there exists a real number 
p = p(t) < 1 such that given $ and F in gL(^B{0,r), then ||m($)-m(F)|| < p||<I>-F||. 



By Lemma (4.8), there exists ^ and C such that rj = ^ + e ■ t(<l'(C)(C) ^ ^c), and 
rj — ( -\- e ■ i(F(C)(^) — Tc). First we compare the ice boundaries |^ — CI to their 
temperatures |$(0(C) " r(C)(C)l- 

IC-CI<e-t|$(e)(e)-r(c)(C)l 

< 6 . t iimio - mm + mm - + mm - mm 

< e • t||<i> - F|| by the definition 

+ ret\^ - CI since for each ^, F(C) G S and | |F| | < r 
+ Let\^-C\) since Fe^^L 
<e-t||<i>-F||+e(L + r)|e-CI 
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Solving for |^ — CI '^6 gst the inequality: 



Let 



and define: 



ie-CI<^^^ll<J>-r|| (4.6) 



e < o.r. , .^ (4-7) 



Then 



2{Lr + L + r) 



Si<52<-. (4.10) 



2 

We now estimate the graph transform map m: 
|m($)(r,)(y)-m(r)(r,)(y)| 

= moiy) + t-[Q- ^iy) ■ (1 - -iB + c)m{y) + cm - A] 



miy) -t-[Q- s{y) ■ (1 - aiOiy)) -{B + C)m{y) + Cm - A] 
miy) + t-[Q- s{y) ■ (1 - a(C)(y)) -{B + C)m{y) + cm - A] 



r(C)(y) -t-[Q- s{y) ■ (1 - aiOiy)) -{B + C)T{Q){y) + CT(C) - A]\ 



Since t < and therefore l — t{B + C) > 0, then the estimate above continues 

as 

< [1 - tiB + c)]\miy) ~ m{y)\ + 1 ■ cm) - m\ 
+ [1 - t(B + c)]\m{y) - miy)\ + 1 ■ c\m - rxoi 

+ \Q-s{y)-{a{0{y)'O.{C){y))\ 



Using inequality (4.6) we estimate the third through the fifth terms of the above 
inequality: 

< [l~t{B + C)+tC]\\<^-T\\ 

+ [1 - t{B + C)]{t5i)\\^ - r|| + {tC){t5i)\\^ ~ T\\ 

+ t52\\^-V\\ 

<{l-tB + t[{l - tB)5i + <52]) 11$ - r|| 
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By the choice of e m (4.7) above, we get the inequality (4.10), and so the inequality 
above continues as: 

< {{l-tB)+t{S2 + S2))\\^-T\\ 

< 1 • ll$-rii 



Therefore, for any fixed t, < t < -g^, if p := ((1 - tB) + t{S2 + S2)), then 
< p < 1, and we showed that: 

||m($) - m(r)|| < p||4> - r||. 

That is, the map to is a contraction. □ 



We finish the proof of the existence of invariant manifold (4.3): 

Proof. For e as in the previous lemma, we showed that 

m:gLnB{0,r)^gcnB{0,r) 

is a contraction in a closed set of a Banach space. Therefore, there exists a unique 
fixed point $* such that to($*) = <&*.□ 

Corollary 4.11. The invariant manifold $* is within 0(e) of the equilibrium 
set T* 

Proof. Since to($*) = then the following holds: 
where £, = rj + fc$. (77). 

Using the facts that F{T*{S,),C} — ^nd the reverse triangle inequality we find: 



> t{B + C) ■ Wii) - T*mB - tB\\W(^) -T^)\\ts 

>tB\\^*{o~T*mB 



In the last step of the above estimate, we used the fact that 
Therefore, we arrive at the following estimate: 



< 



1 



<^-\^*{r, + k^'{r])){Tl + k^'{ri))~T,\ 
< 



B 
er 

B 
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recall that here, r is the bound on the temperature profiles. □ 
Corollary 4.12. For e small, the ice free earth is unstable. 
Proof. Since ||$* - T*|| < then \<^*iv){v) " T*{ri){r])\ < ^ as weh. So that 

if e < f [Tc - r*(l)(l)], then <E>*(1)(1) < T^, and so, ice will form at the pole and 

advances toward the small ice cap equilibrium. See figure (??) for the graph of the 

iceline local equilibrium temperature T* (r/) (ry) . 

Therefore, if r/(0) = 1, as t — > oo, the ice line, advances equatorward, and 
the planet evolves toward a non ice free earth. □ 

5. Concluding Discussion. 

5.1. Future Directions. There are several directions that one can explore based 
on this model. One immediate improvement is to compute the invariant manifold ex- 
plicitly as done by Foias, Sell and Titi [21]. Another immediate improvement on this 
model is an extension to the southern hemisphere with two ice line and a non sym- 
metric transport. Another direction is to explore how the change in the greenhouse 
gas components, which is the term A + BT{y), affects the radiative forcings. A work 
by Andrew Hogg |T5] relates the evolution of temperature with that of carbondioxide 
as a response to the solar input variations caused by the Milankovitch cycle. We 
are interested in the possibility of coupling Budyko's model with the Hogg's model 
to understand the glacial cycles in the quaternary period. While North explores a 
similar model with only a diffusion transport [5], [3], [3], the model discussed in this 
paper could be improved by including some difussion and averaging in the transport 
term. Such inclusion necesitates the consideration for the planet's heat capacity and 
a further explanation of the parameter e. 

5.2. Conclusion. We have shown in this paper the existence of a center stable 
manifold in an energy balance model with ice albedo feedback, featuring a dynamic 
iceline. The existence of such invariant manifold explains the numerical experiments 
presented as animations and allows for qualitative analysis of the small icecap stability. 
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